# P` L ` O ` T ` T ` I ` N ` G

# we are using this for plot: ###### XXXX_-XX_010





#######		#######		#######		#######		#######
#######		#######		#######		#######		#######
#######		#######		#######		#######		#######
#######		#######		#######		#######		#######
#######		#######		#######		B I A S
#######		#######		#######		B I A S
#######		#######		#######		B I A S
#######		#######		#######		B I A S
#######		#######		#######		#######		#######
#######		#######		#######		#######		#######
#######		#######		#######		#######		#######
#######		#######		#######		#######		#######
#######		#######		#######		#######		#######



rm(list=ls())
setwd("/Users/lleemann/Desktop/Replication Files/Simulation/Resulting data")
load("Data_Sample_1000.RData")

objects()

par(family="CMU Serif", mfrow=c(1,2))
#par(family="CMU Serif")
plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,6,2])-1), xlim=c(-1,1), ylim=c(-0.05,0.35), pch=16, col="grey61", cex=2.6, ylab="Bias", xlab="Correlation of Errors", main=expression(paste("Estimate of ", beta[21])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,6,2])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,6,2])-1),pch=17, col="black", cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,6,6])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,6,6])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,6,6])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,6,10])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,6,10])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,6,10])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,6,14])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,6,14])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,6,14])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,6,18])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,6,18])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,6,18])-1),pch=17, col="black", cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,6,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,6,6])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,6,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,6,10])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,6,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,6,14])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,6,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,6,18])-1), lty=2, col="grey61", lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,6,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,6,6])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,6,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,6,10])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,6,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,6,14])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,6,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,6,18])-1), lty=2, col="grey80", lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,6,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,6,6])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,6,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,6,10])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,6,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,6,14])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,6,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,6,18])-1), lty=2, col="black", lwd=2)


legend(.4,.3, legend=c("Strat", "Heck", "SSa"), col=c("grey61","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")


### 2nd plot

plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,5,2])-1), xlim=c(-1,1), ylim=c(-0.05,0.35), pch=16, col="grey61", cex=2.6, ylab="Bias", xlab="Correlation of Errors", main=expression(paste("Estimate of ", beta[22])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,5,2])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,5,2])-1),pch=17, col="black", cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,5,6])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,5,6])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,5,6])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,5,10])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,5,10])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,5,10])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,5,14])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,5,14])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,5,14])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,5,18])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,5,18])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,5,18])-1),pch=17, col="black", cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,5,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,5,6])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,5,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,5,10])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,5,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,5,14])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,5,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,5,18])-1), lty=2, col="grey61", lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,5,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,5,6])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,5,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,5,10])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,5,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,5,14])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,5,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,5,18])-1), lty=2, col="grey80", lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,5,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,5,6])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,5,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,5,10])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,5,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,5,14])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,5,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,5,18])-1), lty=2, col="black", lwd=2)



legend(.4,.3, legend=c("Strat", "Heck", "SSa"), col=c("grey61","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")

####

#### 


####

#### 


####

#par(family="CMU Serif", mfrow=c(1,2))
par(family="CMU Serif")
plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,3,2])-1), xlim=c(-1,1), ylim=c(-0.65,0.25), pch=16, col="grey61", cex=2.6, ylab="Bias", xlab="Correlation of Errors", main=expression(paste("Estimate of ", beta[12])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,3,2])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,3,2])-1),pch=17, col="black", cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,3,6])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,3,6])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,3,6])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,3,10])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,3,10])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,3,10])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,3,14])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,3,14])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,3,14])-1),pch=17, col="black", cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,3,18])-1), pch=16, col="grey61", cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,3,18])-1),pch=15, col="grey80", cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,3,18])-1),pch=17, col="black", cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,3,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,3,6])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,3,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,3,10])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,3,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,3,14])-1), lty=2, col="grey61", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,3,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,3,18])-1), lty=2, col="grey61", lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,3,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,3,6])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,3,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,3,10])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,3,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,3,14])-1), lty=2, col="grey80", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,3,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,3,18])-1), lty=2, col="grey80", lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,3,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,3,6])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,3,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,3,10])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,3,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,3,14])-1), lty=2, col="black", lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,3,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,3,18])-1), lty=2, col="black", lwd=2)


legend(.4,.23, legend=c("Strat", "Heck", "SSa"), col=c("grey61","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")


